clear; clc; close all;

alpha = linspace(-6, 6, 1200);
eta = linspace(-6, 6, 1200);
[Alpha, Eta] = meshgrid(alpha, eta);

N = numel(Alpha);
ClassMapVec = zeros(N,1); % 线性存储分类结果

parfor idx = 1:N
    a = Alpha(idx);
    e = Eta(idx);
    
    m1 = a*exp(-e) - 0.5;
    m2 = -0.5*a*exp(-e) + 4.5;
    m3 = 8*a*exp(-e) - 8;
    
    roots_nonzero = roots([1 m1 m2 m3]);
    
    nPos = sum(real(roots_nonzero) > 1e-8);
    nNeg = sum(real(roots_nonzero) < -1e-8);
    
    if nPos == 3 && nNeg == 0
        ClassMapVec(idx) = 1; % 绿 3P0N
    elseif nPos == 2 && nNeg == 1
        ClassMapVec(idx) = 2; % 红 2P1N
    elseif nPos == 1 && nNeg == 2
        ClassMapVec(idx) = 3; % 黄 1P2N
    elseif nPos == 0 && nNeg == 3
        ClassMapVec(idx) = 4; % 蓝 0P3N（稳定）
    else
        ClassMapVec(idx) = 0; % 其他
    end
end

ClassMap = reshape(ClassMapVec, size(Alpha));

% 绘图
figure;
imagesc(alpha, eta, ClassMap);
set(gca,'YDir','normal');
colormap([0 1 0; 1 0 0; 1 1 0; 0 0 1]); % 绿, 红, 黄, 蓝
colorbar('Ticks',1:4, ...
    'TickLabels',{'3P0N','2P1N','1P2N','0P3N'});
caxis([0.5, 4.5]);
xlabel('\alpha', 'FontSize', 12);
ylabel('\eta', 'FontSize', 12);

function dx = mCNN(x, alpha)
    A = [-alpha, 0, 0;
          0,     0, 0;
          0,     0, 0];
    B = [1, -4, -3.5;
         0,  1,  2;
        -1, -4, 1.5];
    X = x(1:3);
    Q = x(4:6);
    M = exp(-Q);
    Y = M .* X;
    dXdt = -X + A*Y + B*X;
    dQdt = X;
    dx = [dXdt; dQdt];
end